feedback
Losing track of variable definitions
- \(\alpha\) and \(\beta\) are confusing. Need to think about how to frame them better in the slides.
- Options include:
- footnotes/box in the corner with definitions for reminder
- always writing them as the ratios
- I wonder if it would be clearer to abandom \(S_1, S_0\) and instead spell out \(\Pr(\text{Symptomatic}\,|\,\text{untested})\), \(\Pr(\text{Asymptomatic}\,|\,\text{untested})\) (though clarify that’s not exactly what they mean)
Evan’s question on dependency of \(\alpha\) and \(\beta\)
- \(\alpha = \dfrac{P(\text{test}_+| \text{untested}, S_1)}{P(\text{test}_+|\text{tested})}\)
- \(\beta = \dfrac{P(\text{test}_+| \text{untested}, S_0)}{P(\text{test}_+|\text{tested})}\)
- denominator is not treated as a random variable – this is the test positivity we observe
Notes on Independence (Pre-melding)
The number of events included here makes thinking about dependencies difficult, particularly since \(\alpha\) and \(\beta\) are related, but we sample them independently. It’s more clear to ignore the denominator since \(P(\text{test}_+|\text{tested})\) is just the observed test positivity and is a constant.
Observing \(P(\text{test}_+|\text{untested},S_0)\) doesn’t tell us anything directly about \(P(\text{test}_+|\text{untested},S_1)\). For example, if \(P(\text{test}_+|\text{untested},S_0)=1\), \(P(\text{test}_+|\text{untested},S_1)\) could be zero, or it could be 1, or anything in between.
If we had information about \(P(\text{test}_+|\text{untested}), P(S_1),\) and \(P(S_0)\), we could relate these quantities as
\[P(\text{test}_+|\text{untested}) = P(\text{test}_+, S_0|\text{untested}) P(S_0) + P(\text{test}_+, S_0|\text{untested})P(S_1).\]
We can also see their relationship with \(P(S_1|\text{untested}) = 1 - P(S_0|\text{untested})\) by considering
\[P(\text{test}_+|\text{untested}, S_0) = \dfrac{P(\text{test}_+, \text{untested}, S_0)}{P(\text{untested},S_0)} = \dfrac{P(\text{test}_+, S_0| \text{untested})}{P(S_0|\text{untested})},\] and similarly
\[P(\text{test}_+|\text{untested}, S_1) = \dfrac{P(\text{test}_+, S_1| \text{untested})}{P(S_1|\text{untested})} = \dfrac{P(\text{test}_+, S_1| \text{untested})}{1-P(S_0|\text{untested})},\]
but that said, with information on \(P(S_1|\text{untested})\) alone, we will would not be able to relate \(P(\text{test}_+|\text{untested}, S_0)\) and \(P(\text{test}_+|\text{untested}, S_1)\) directly because we wouldn’t know how the numerators \(P(\text{test}_+, S_0| \text{untested}), \;\; P(\text{test}_+, S_1| \text{untested})\) relate.
\[M(\theta) = \dfrac{\beta(1- P(S_1|\text{untested}))}{\beta(1-P(S_1|\text{untested})) + \alpha P(S_1|\text{untested})}\]
Summary: \(\alpha\) and \(\beta\) are dependent, but they are dependent on information we don’t have, so without this knowledge we sample them independently.
Notes on Independence (Post-melding)
During melding, we resample from \(\alpha, \beta, P(S_1|\text{untested})\) with the same weights, placing more weight on observations that yield observations with high density in the pooled distribution. If we compute \(M(\theta)\) with the constrained distributions of \(\theta\), the resulting distribution is the melded distribution.
When we sample from them from the melded priors for use in the probabilistic bias analysis, we account for this dependency by sampling them together.
source(here('analysis/base_functions/get_melded.R'))
source(here('analysis/base_functions/base_functions.R'))
prior_params <- list(
alpha_mean = .95,
alpha_sd = 0.08,
alpha_bounds = NA,
# alpha_bounds = c(.8,1),
beta_mean = .15,
beta_sd =.09,
beta_bounds = NA,
# beta_bounds = c(0.002, 0.4),
s_untested_mean = .03,
s_untested_sd = .0225,
# s_untested_bounds = c(0.0018, Inf),
s_untested_bounds = NA,
p_s0_pos_mean = .4,
# p_s0_pos_sd = .1225,
p_s0_pos_sd = .1,
p_s0_pos_bounds = NA,
# p_s0_pos_bounds = c(.25, .7),
pre_nsamp = 1e6,
post_nsamp = 1e5)
melded <- do.call(get_melded, prior_params)
melded$pre_melding %>%
slice_sample(n=1e5) %>%
mutate(source = "Before Melding") %>%
bind_rows(melded$post_melding) %>%
mutate(source = ifelse(is.na(source), "After Melding", source)) %>%
ggplot(aes(x=alpha,y=beta,color = source)) +
geom_point(alpha =.2,size=.5)+
scale_color_manual(values = c(
"Before Melding"= "#5670BF",
"Induced" = "#B28542",
"After Melding" = "#418F6A")) +
labs(color = "",
x = TeX("$\\alpha$"),
y = TeX("$\\beta$"),
title = TeX("Relationship between $\\alpha$ and $\\beta$ Before and After Melding")) +
theme_c(legend.position="top") +
guides(color = guide_legend(override.aes = list(size = 5, alpha = 1)))melded$post_melding %>%
ggplot(aes(x=alpha,y=beta)) +
geom_point(alpha =.2,size=.5)
melded <- do.call(get_melded, prior_params)
post <- melded$post_melding %>%
select(P_A_testpos) %>%
pivot_longer(everything())
melded$post_melding %>%
mutate(induced_new = est_P_A_testpos(P_S_untested, alpha, beta)) %>%
select(induced_new) %>%
pivot_longer(everything()) %>%
bind_rows(post) %>%
ggplot(aes(x=value, fill =name)) +
geom_density(alpha=.9)
melded$post_melding %>%
mutate(induced_new = est_P_A_testpos(P_S_untested, alpha, beta)) %>%
filter(induced_new != P_A_testpos)
melded_long <- reformat_melded(melded$post_melding,
melded$pre_melding,
pre_nsamp = prior_params$pre_nsamp,
p_s0_pos_mean = prior_params$p_s0_pos_mean,
p_s0_pos_sd = prior_params$p_s0_pos_sd,
p_s0_pos_bounds = NA)
################################
# FIRST WITHOUT MELDING
################################
melded_long %>%
filter(type == "Before Melding") %>%
ggplot(aes(x=value, fill = type)) +
geom_density( alpha = .8, show.legend=FALSE) +
facet_wrap(~name, labeller=as_labeller(TeX, default = label_parsed), ncol=4) +
theme_c(legend.position = "top") +
scale_fill_manual(values = list("Before Melding"= "#5670BF", "Induced" = "#B28542")) +
labs(fill = "")
plt2 <- melded_long %>%
filter(type!="After Melding" & name == "$P(S_0|test+,untested)$") %>%
ggplot(aes(x=value, fill = type)) +
geom_density( alpha = .8, show.legend=TRUE) +
facet_wrap(~name, labeller=as_labeller(TeX, default = label_parsed), ncol=4) +
theme_c(legend.position = "right",
legend.text = element_text(size = 16)) +
scale_fill_manual(values = list("Before Melding"= "#5670BF", "Induced" = "#B28542")) +
labs(fill = "")
cowplot::plot_grid(plt1,plt2, ncol=1, rel_widths = c(1,.9))
ggsave(width=9, height = 8, "./presentation/figure/plot_no_melding.jpeg", dpi=300)
################################
# NOW WITH MELDING
################################
plt1 <- melded_long %>%
filter(name != "$P(S_0|test+,untested)$") %>%
ggplot(aes(x=value, fill = type)) +
geom_density( alpha = .8, show.legend=FALSE) +
facet_wrap(~name, labeller=as_labeller(TeX, default = label_parsed), ncol=4) +
theme_c(legend.position = "top") +
scale_fill_manual(values = c("Before Melding"= "#5670BF", "Induced" = "#B28542",
"After Melding" = "#418F6A")) +
labs(fill = "")
plt2 <- melded_long %>%
filter( name == "$P(S_0|test+,untested)$") %>%
ggplot(aes(x=value, fill = type)) +
geom_density( alpha = .8, show.legend=TRUE) +
facet_wrap(~name, labeller=as_labeller(TeX, default = label_parsed), ncol=4) +
theme_c(legend.position = "right",
legend.text = element_text(size = 16)) +
scale_fill_manual(values = c("Before Melding"= "#5670BF", "Induced" = "#B28542",
"After Melding" = "#418F6A")) +
labs(fill = "")
cowplot::plot_grid(plt1,plt2, ncol=1, rel_widths = c(1,.9))Evan’s question on why \(\alpha\) doesn’t change pre/post melding
Intuition: we form weights to sample from the log-pooled prior of \(\phi\), so since these weights aren’t heavily dependent on \(\alpha\), we are basically just randomly sampling from \(\alpha\).
prior_params <- list(
alpha_mean = .95,
alpha_sd = 0.08,
alpha_bounds = NA,
# alpha_bounds = c(.8,1),
beta_mean = .15,
beta_sd =.09,
beta_bounds = NA,
# beta_bounds = c(0.002, 0.4),
s_untested_mean = .03,
s_untested_sd = .0225,
# s_untested_bounds = c(0.0018, Inf),
s_untested_bounds = NA,
p_s0_pos_mean = .4,
# p_s0_pos_sd = .1225,
p_s0_pos_sd = .1,
p_s0_pos_bounds = NA,
# p_s0_pos_bounds = c(.25, .7),
pre_nsamp = 1e6,
post_nsamp = 1e5)
get_theta <- function(alpha_mean = 0.9,
alpha_sd = 0.04,
alpha_bounds = NA,
beta_mean = .15,
beta_sd =.09,
beta_bounds = NA,
s_untested_mean = .025,
s_untested_sd = .0225,
s_untested_bounds = NA,
p_s0_pos_mean = .4,
p_s0_pos_sd = .1225,
p_s0_pos_bounds = NA,
pre_nsamp = 1e4,
post_nsamp = 1e3,
include_corrected = TRUE,
bde = FALSE) {
tibble(alpha = sample_gamma_density(pre_nsamp,
mean = alpha_mean,
sd = alpha_sd,
bounds = alpha_bounds),
beta= sample_beta_density(pre_nsamp,
mean = beta_mean,
sd = beta_sd,
bounds = beta_bounds),
P_S_untested = sample_beta_density(pre_nsamp,
mean = s_untested_mean,
sd = s_untested_sd,
bounds = s_untested_bounds)) %>%
mutate(phi_induced = est_P_A_testpos(P_S_untested = P_S_untested,
alpha = alpha,
beta=beta))
}
tibble(alpha =seq(.7,1.3, length =100),
beta = .15,
P_S_untested = .03,
induced = (beta * (1 - P_S_untested)) / (( beta * (1 - P_S_untested)) + (alpha * P_S_untested))) %>%
ggplot(aes(x=alpha, y = induced)) +
geom_line() +
theme_c() +
labs(x= TeX("$\\alpha$"),
y= TeX("$M(\\theta)$"),
title = TeX("$M(\\theta)$ Across Values of $\\alpha$"),
subtitle = TeX("$P(S_1|untested),\\beta$ are Fixed"))+
ylim(.4, .95)
tibble(beta =seq(.05,.45, length =100),
alpha = .95,
P_S_untested = .03,
induced = (beta * (1 - P_S_untested)) / (( beta * (1 - P_S_untested)) + (alpha * P_S_untested))) %>%
ggplot(aes(x=beta, y = induced)) +
geom_line()+
theme_c() +
labs(x= TeX("$\\beta$"),
y= TeX("$M(\\theta)$"),
title = TeX("$M(\\theta)$ Across Values of $\\beta$"),
subtitle = TeX("$P(S_1|untested),\\alpha$ are Fixed")) +
ylim(.4, .95)
tibble(P_S_untested =seq(.05,.18, length =100),
alpha = .95,
beta=.15,
induced = (beta * (1 - P_S_untested)) / (( beta * (1 - P_S_untested)) + (alpha * P_S_untested))) %>%
ggplot(aes(x=P_S_untested, y = induced)) +
geom_line()+
theme_c() +
labs(x= TeX("$\\beta$"),
y= TeX("$M(\\theta)$"),
title = TeX("$M(\\theta)$ Across Values of $P(S_1|untested)$"),
subtitle = TeX("$\\beta,\\alpha$ are Fixed")) +
ylim(.4, .95)tibble(alpha =seq(.7,1.3, length =100),
beta = .15,
P_S_untested = .03,
induced = (beta * (1 - P_S_untested)) / (( beta * (1 - P_S_untested)) + (alpha * P_S_untested))) %>%
ggplot(aes(x=alpha, y = induced)) +
geom_line() +
theme_c() +
labs(x= TeX("$\\alpha$"),
y= TeX("$M(\\theta)$"),
title = TeX("$M(\\theta)$ Across Values of $\\alpha$"),
subtitle = TeX("$P(S_1|untested),\\beta$ are Fixed"))+
ylim(.4, .95)In Figure 1, we see increasing the mean of alpha corresponds to less density toward 1 and greater density at lower values of \(\phi\), but the effect is relatively small.
res <- map_df(seq(.8,1.3, by = .1), ~ {
params <- prior_params
params$alpha_mean <- .x
theta <- do.call(get_theta, params)
theta %>%
mutate(alpha_mean = .x)
} )
res %>%
mutate(alpha_mean = factor(alpha_mean)) %>%
ggplot(aes(x=phi_induced, fill = alpha_mean)) +
geom_density(alpha =.5) +
viridis::scale_fill_viridis(option = "mako", discrete=TRUE, begin=.1, end = .9,direction=-1) +
theme_c(legend.position="right",
legend.text = element_text(size = 10),
legend.title = element_text(size = 14)) +
labs(title = TeX("Induced Distribution on $\\phi$ by Mean of Alpha"),
x = TeX("$\\phi^{induced}$"),
fill = TeX("Mean of $\\alpha$")) +
guides(fill = guide_legend(override.aes = list(alpha =1)))Figure 1:
In particular, the effect of \(\alpha\) is limited by the small values of \(P(S_1|\text{untested})\). When \(P(S_1|untested)\) increases, \(\alpha\) makes a larger difference in the value of \(M(\theta)\) ( 2).
Red line corresponds to where \(P(S_1|\text{untested})\) is closest to its mean.
res <- map_df(seq(0.001, .2, length = 50), ~ {
tibble(alpha =seq(.7,1.3, length =100),
beta = .15,
P_S_untested = .x,
induced = (beta * (1 - P_S_untested)) / (( beta * (1 - P_S_untested)) + (alpha * P_S_untested)))
})
realistic <- res %>%
filter(abs(P_S_untested - 0.02536735) < 1e-4)
res %>%
ggplot(aes(x=alpha, y = induced, color = P_S_untested, group =P_S_untested )) +
geom_line() +
geom_line(data = realistic,
aes(x=alpha,
y = induced,
color = P_S_untested,
group =P_S_untested ),
color = 'red',
linewidth =1.5) +
labs(x= TeX("$\\alpha$"),
y= TeX("$M(\\theta)$"),
title = TeX("$M(\\theta)$ Across Values of $\\alpha$"),
subtitle = TeX("$P(S_1|untested),\\beta$ are Fixed"),
color = TeX("$P(S_1|untested)$"))+
viridis::scale_color_viridis(option="magma", end=.9, direction = -1)+
theme_c(legend.position="right",
legend.text = element_text(size = 10),
legend.title = element_text(size = 14))
res <- map_df(seq(0.001, .2, length = 50), ~ {
tibble(alpha =.95,
beta = seq(.03,.4, length =100),
P_S_untested = .x,
induced = (beta * (1 - P_S_untested)) / (( beta * (1 - P_S_untested)) + (alpha * P_S_untested)))
})
realistic <- res %>%
filter(abs(P_S_untested - 0.02536735) < 1e-4)
res %>%
ggplot(aes(x=beta, y = induced, color = P_S_untested, group =P_S_untested )) +
geom_line() +
geom_line(data = realistic,
aes(x=beta,
y = induced,
color = P_S_untested,
group =P_S_untested ),
color = 'red',
linewidth =1.5) +
labs(x= TeX("$\\beta$"),
y= TeX("$M(\\theta)$"),
title = TeX("$M(\\theta)$ Across Values of $\\beta$"),
subtitle = TeX("$P(S_1|untested),\\alpha$ are Fixed"),
color = TeX("$P(S_1|untested)$"))+
viridis::scale_color_viridis(option="magma", end=.9, direction = -1)+
theme_c(legend.position="right",
legend.text = element_text(size = 10),
legend.title = element_text(size = 14))Figure 2:
Comparing Melded Distributions by Mean of \(\alpha\)
prior_params <- list(
alpha_mean = .95,
alpha_sd = 0.08,
alpha_bounds = NA,
# alpha_bounds = c(.8,1),
beta_mean = .15,
beta_sd =.09,
beta_bounds = NA,
# beta_bounds = c(0.002, 0.4),
s_untested_mean = .03,
s_untested_sd = .0225,
# s_untested_bounds = c(0.0018, Inf),
s_untested_bounds = NA,
p_s0_pos_mean = .4,
# p_s0_pos_sd = .1225,
p_s0_pos_sd = .1,
p_s0_pos_bounds = NA,
# p_s0_pos_bounds = c(.25, .7),
pre_nsamp = 1e5,
post_nsamp = 1e4)
res <- map_df(seq(.5, 2, length =5),
~ {
params <- prior_params
params$alpha_mean <- .x
melded <- do.call(get_melded, params)
reformat_melded(melded$post_melding,
melded$pre_melding,
params$pre_nsamp,
p_s0_pos_mean = params$p_s0_pos_mean,
p_s0_pos_sd = params$p_s0_pos_sd,
p_s0_pos_bounds = params$p_s0_pos_bounds) %>%
mutate(alpha_mean = .x)
})
res %>%
mutate(alpha_mean = factor(round(alpha_mean,1))) %>%
# filter(name == "$\\alpha$") %>%
ggplot(aes(x = value, fill = type)) +
geom_density(alpha=.8) +
facet_grid(name~alpha_mean,
scales = "free_y",
labeller=as_labeller(TeX,
default=label_parsed)) +
scale_fill_manual(values = c("Before Melding"= "#5670BF", "Induced" = "#B28542",
"After Melding" = "#418F6A")) +
theme_c() +
labs(fill="")Comparing Melded Distributions by Mean of \(\alpha\) with Increased \(\alpha\) SD
Increased SD
res <- map_df(seq(.5, 2, length =5),
~ {
params <- prior_params
params$alpha_mean <- .x
params$alpha_sd <- .3
# params$s_untested_mean <- .2
melded <- do.call(get_melded, params)
reformat_melded(melded$post_melding,
melded$pre_melding,
params$pre_nsamp,
p_s0_pos_mean = params$p_s0_pos_mean,
p_s0_pos_sd = params$p_s0_pos_sd,
p_s0_pos_bounds = params$p_s0_pos_bounds) %>%
mutate(alpha_mean = .x)
})
res %>%
mutate(alpha_mean = factor(round(alpha_mean,1))) %>%
# filter(name == "$\\alpha$") %>%
ggplot(aes(x = value, fill = type)) +
geom_density(alpha=.8) +
facet_grid(name~alpha_mean,
scales = "free_y",
labeller=as_labeller(TeX,
default=label_parsed)) +
scale_fill_manual(values = c("Before Melding"= "#5670BF", "Induced" = "#B28542",
"After Melding" = "#418F6A")) +
theme_c() +
labs(fill="")Decreased SD
res <- map_df(seq(.5, 2, length =5),
~ {
params <- prior_params
params$alpha_mean <- .x
params$alpha_sd <- .04
# params$s_untested_mean <- .2
melded <- do.call(get_melded, params)
reformat_melded(melded$post_melding,
melded$pre_melding,
params$pre_nsamp,
p_s0_pos_mean = params$p_s0_pos_mean,
p_s0_pos_sd = params$p_s0_pos_sd,
p_s0_pos_bounds = params$p_s0_pos_bounds) %>%
mutate(alpha_mean = .x)
})
res %>%
mutate(alpha_mean = factor(round(alpha_mean,1))) %>%
# filter(name == "$\\alpha$") %>%
ggplot(aes(x = value, fill = type)) +
geom_density(alpha=.8) +
facet_grid(name~alpha_mean,
scales = "free_y",
labeller=as_labeller(TeX,
default=label_parsed)) +
scale_fill_manual(values = c("Before Melding"= "#5670BF", "Induced" = "#B28542",
"After Melding" = "#418F6A")) +
theme_c() +
labs(fill="")Comparing Melded Distributions by Mean of \(\alpha\) with Increased \(P(S_1|\text{untested})\) Mean
Upping \(P(S_1|\text{untested})\) to be .5
Unreasonable, but just to see
res <- map_df(seq(.5, 2, length =5),
~ {
params <- prior_params
params$alpha_mean <- .x
params$s_untested_mean <- .3
# params$s_untested_mean <- .2
melded <- do.call(get_melded, params)
reformat_melded(melded$post_melding,
melded$pre_melding,
params$pre_nsamp,
p_s0_pos_mean = params$p_s0_pos_mean,
p_s0_pos_sd = params$p_s0_pos_sd,
p_s0_pos_bounds = params$p_s0_pos_bounds) %>%
mutate(alpha_mean = .x)
})
res %>%
mutate(alpha_mean = factor(round(alpha_mean,1))) %>%
# filter(name == "$\\alpha$") %>%
ggplot(aes(x = value, fill = type)) +
geom_density(alpha=.8) +
facet_grid(name~alpha_mean,
scales = "free_y",
labeller=as_labeller(TeX,
default=label_parsed)) +
scale_fill_manual(values = c("Before Melding"= "#5670BF", "Induced" = "#B28542",
"After Melding" = "#418F6A")) +
theme_c() +
labs(fill="")Upping \(P(S_1|\text{untested})\) to be .5
res <- map_df(seq(.5, 2, length =5),
~ {
params <- prior_params
params$alpha_mean <- .x
params$s_untested_mean <- .5
# params$s_untested_mean <- .2
melded <- do.call(get_melded, params)
reformat_melded(melded$post_melding,
melded$pre_melding,
params$pre_nsamp,
p_s0_pos_mean = params$p_s0_pos_mean,
p_s0_pos_sd = params$p_s0_pos_sd,
p_s0_pos_bounds = params$p_s0_pos_bounds) %>%
mutate(alpha_mean = .x)
})
res %>%
mutate(alpha_mean = factor(round(alpha_mean,1))) %>%
# filter(name == "$\\alpha$") %>%
ggplot(aes(x = value, fill = type)) +
geom_density(alpha=.8) +
facet_grid(name~alpha_mean,
scales = "free_y",
labeller=as_labeller(TeX,
default=label_parsed)) +
scale_fill_manual(values = c("Before Melding"= "#5670BF", "Induced" = "#B28542",
"After Melding" = "#418F6A")) +
theme_c() +
labs(fill="")Comparing Melded Distributions by Mean of \(\alpha\) with Increased \(\alpha\) SD
Increased to .15
prior_params <- list(
alpha_mean = .95,
alpha_sd = 0.08,
# alpha_sd = 0.15,
alpha_bounds = NA,
# alpha_bounds = c(.8,1),
beta_mean = .15,
beta_sd =.09,
beta_bounds = NA,
# beta_bounds = c(0.002, 0.4),
s_untested_mean = .03,
s_untested_sd = .0225,
# s_untested_bounds = c(0.0018, Inf),
s_untested_bounds = NA,
p_s0_pos_mean = .4,
# p_s0_pos_sd = .1225,
p_s0_pos_sd = .1,
p_s0_pos_bounds = NA,
# p_s0_pos_bounds = c(.25, .7),
pre_nsamp = 1e5,
post_nsamp = 1e4)
res <- map_df(seq(.5, 2, length =5),
~ {
params <- prior_params
params$alpha_mean <- .x
params$alpha_sd <- .15
params$s_untested_mean <- .03
melded <- do.call(get_melded, params)
# melded$post_melding %>%
# mutate(alpha_mean = .x) %>%
# cbind(alpha_before = melded$pre_melding$alpha)
reformat_melded(melded$post_melding,
melded$pre_melding,
params$pre_nsamp,
p_s0_pos_mean = params$p_s0_pos_mean,
p_s0_pos_sd = params$p_s0_pos_sd,
p_s0_pos_bounds = params$p_s0_pos_bounds) %>%
mutate(alpha_mean = .x)
})
res %>%
mutate(alpha_mean = factor(round(alpha_mean,1))) %>%
# filter(name == "$\\alpha$") %>%
ggplot(aes(x = value, fill = type)) +
geom_density(alpha=.8) +
facet_grid(name~alpha_mean,
scales = "free_y",
labeller=as_labeller(TeX,
default=label_parsed)) +
scale_fill_manual(values = c("Before Melding"= "#5670BF", "Induced" = "#B28542",
"After Melding" = "#418F6A")) +
theme_c() +
labs(fill="")#
# res <- map_df(seq(.5, 2, length =5),
# ~ {
# params <- prior_params
# params$alpha_mean <- .x
# params$alpha_sd <- .15
#
# params$s_untested_mean <- .2
# melded <- do.call(get_melded, params)
# # melded$post_melding %>%
# # mutate(alpha_mean = .x) %>%
# # cbind(alpha_before = melded$pre_melding$alpha)
#
# reformat_melded(melded$post_melding,
# melded$pre_melding,
# params$pre_nsamp,
# p_s0_pos_mean = params$p_s0_pos_mean,
# p_s0_pos_sd = params$p_s0_pos_sd,
# p_s0_pos_bounds = params$p_s0_pos_bounds) %>%
# mutate(alpha_mean = .x)
#
# })
#
#
# res %>%
# mutate(alpha_mean = factor(round(alpha_mean,1))) %>%
# # filter(name == "$\\alpha$") %>%
# ggplot(aes(x = value, fill = type)) +
# geom_density(alpha=.8) +
# facet_grid(name~alpha_mean,
# scales = "free_y",
# labeller=as_labeller(TeX,
# default=label_parsed)) +
# scale_fill_manual(values = c("Before Melding"= "#5670BF", "Induced" = "#B28542",
# "After Melding" = "#418F6A")) +
# theme_c() +
# labs(fill="")
# Increased to .25
res <- map_df(seq(.5, 2, length =5),
~ {
params <- prior_params
params$alpha_mean <- .x
params$alpha_sd <- .25
params$s_untested_mean <- .03
melded <- do.call(get_melded, params)
# melded$post_melding %>%
# mutate(alpha_mean = .x) %>%
# cbind(alpha_before = melded$pre_melding$alpha)
reformat_melded(melded$post_melding,
melded$pre_melding,
params$pre_nsamp,
p_s0_pos_mean = params$p_s0_pos_mean,
p_s0_pos_sd = params$p_s0_pos_sd,
p_s0_pos_bounds = params$p_s0_pos_bounds) %>%
mutate(alpha_mean = .x)
})
res %>%
mutate(alpha_mean = factor(round(alpha_mean,1))) %>%
# filter(name == "$\\alpha$") %>%
ggplot(aes(x = value, fill = type)) +
geom_density(alpha=.8) +
facet_grid(name~alpha_mean,
scales = "free_y",
labeller=as_labeller(TeX,
default=label_parsed)) +
scale_fill_manual(values = c("Before Melding"= "#5670BF", "Induced" = "#B28542",
"After Melding" = "#418F6A")) +
theme_c() +
labs(fill="")Fixing \(\beta,P(S_1|\text{untested})\)
res <- map_df(seq(.5, 2, length =5),
~ {
params <- prior_params
params$alpha_mean <- .x
params$beta_sd <- 0
params$s_untested_sd <- 0
melded <- do.call(get_melded, params)
# melded$post_melding %>%
# mutate(alpha_mean = .x) %>%
# cbind(alpha_before = melded$pre_melding$alpha)
reformat_melded(melded$post_melding,
melded$pre_melding,
params$pre_nsamp,
p_s0_pos_mean = params$p_s0_pos_mean,
p_s0_pos_sd = params$p_s0_pos_sd,
p_s0_pos_bounds = params$p_s0_pos_bounds) %>%
mutate(alpha_mean = .x)
})
res %>%
mutate(alpha_mean = factor(round(alpha_mean,1))) %>%
# filter(name == "$\\alpha$") %>%
ggplot(aes(x = value, fill = type)) +
geom_density(alpha=.8) +
facet_grid(name~alpha_mean,
scales = "free_y",
labeller=as_labeller(TeX,
default=label_parsed)) +
scale_fill_manual(values = c("Before Melding"= "#5670BF", "Induced" = "#B28542",
"After Melding" = "#418F6A")) +
theme_c() +
labs(fill="")Fixed versus Free Scales
Even allowing for different scales, \(\alpha\) does not change substantially.
Fixed Scale
melded <- do.call(get_melded, prior_params)
post <- melded$post_melding %>%
select(P_A_testpos) %>%
pivot_longer(everything())
melded_long <- reformat_melded(melded$post_melding,
melded$pre_melding,
pre_nsamp = prior_params$pre_nsamp,
p_s0_pos_mean = prior_params$p_s0_pos_mean,
p_s0_pos_sd = prior_params$p_s0_pos_sd,
p_s0_pos_bounds = NA)
################################
# NOW WITH MELDING
################################
plt1 <- melded_long %>%
filter(name != "$P(S_0|test+,untested)$") %>%
ggplot(aes(x=value, fill = type)) +
geom_density( alpha = .8, show.legend=FALSE) +
facet_wrap(~name, labeller=as_labeller(TeX, default = label_parsed), ncol=3) +
theme_c(legend.position = "top") +
scale_fill_manual(values = c("Before Melding"= "#5670BF", "Induced" = "#B28542",
"After Melding" = "#418F6A")) +
labs(fill = "")
plt2 <- melded_long %>%
filter( name == "$P(S_0|test+,untested)$") %>%
ggplot(aes(x=value, fill = type)) +
geom_density( alpha = .8, show.legend=TRUE) +
facet_wrap(~name, labeller=as_labeller(TeX, default = label_parsed), ncol=4) +
theme_c(legend.position = "right",
legend.text = element_text(size = 16)) +
scale_fill_manual(values = c("Before Melding"= "#5670BF", "Induced" = "#B28542",
"After Melding" = "#418F6A")) +
labs(fill = "")
cowplot::plot_grid(plt1,plt2, ncol=1, rel_widths = c(1,.9))Free Scale
################################
# NOW WITH MELDING
################################
plt1 <- melded_long %>%
filter(name != "$P(S_0|test+,untested)$") %>%
ggplot(aes(x=value, fill = type)) +
geom_density( alpha = .8, show.legend=FALSE) +
facet_wrap(~name, labeller=as_labeller(TeX, default = label_parsed), ncol=4, scales="free_y") +
theme_c(legend.position = "top") +
scale_fill_manual(values = c("Before Melding"= "#5670BF", "Induced" = "#B28542",
"After Melding" = "#418F6A")) +
labs(fill = "")
cowplot::plot_grid(plt1,plt2, ncol=1, rel_widths = c(1,.9))Conclusions
The lack of an effect of \(\alpha\) is primarily a result of the functional form of \(M\); because \(\beta\) has a larger effect
M <- function(input1, input2) {
# x <- 0.5
x <- 0.03
input1*(1-x) / (input1*(1-x) + input2*x)
}
M <- Vectorize(M )
inp1 <- seq(0.01,.4, length = 50)
names(inp1) <- inp1
inp2 <- seq(.7,1.3, length = 50)
names(inp2) <- inp2
z <- outer(inp1,inp2, FUN =M)
# 3D plot
par(mar = c(0, 0, 0, 0))
persp(inp1, inp2, z, xlab = "beta", ylab= "alpha", theta = 60, phi = 30)library(plotly)
ax <- list(
title = "alpha"
)
ay <- list(
title = "beta"
)
fig <- plot_ly(z=~z)
fig <- plot_ly( x= ~colnames(z),
y =~ rownames(z),
z=~z)
fig %>% add_surface() %>% layout(scene = list(xaxis = ax,
yaxis = ay))Aaron’s question on notation
Evan’s question on \(P(S_0|\text{test}_+, \text{untested}) \neq P(S_0|\text{test}_+, \text{tested})\)
The fact that \(P(S_0|\text{test}_+, \text{untested}) > P(S_0|\text{test}_+, \text{tested})\) is something important to think about. It has been a while since my initial search for meta-analyses, so another look would be good to get a better sense for the estimates from screening studies and the extent to which they differ.